import numpy as np
import math

# 参数
theta = 120  # 开角（度）
alpha = 1.5  # 坡度（度）
D0 = 110     # 中心深度（m）
L_east_west = 4 * 1852  # 东西宽度（m）
L_north_south = 2 * 1852 # 南北长度（m）

# 覆盖宽度函数（问题2模型）
def coverage_width(D, beta):
    theta_half = math.radians(theta/2)  # 半开角（弧度）
    alpha_rad = math.radians(alpha)
    beta_rad = math.radians(beta)
    # 注意：这里D是当地水深，beta是测线与坡面法向投影夹角（坡面法向？假设坡面法向水平投影为正东方向）
    # 由于西深东浅，坡面法向大致向东（即上坡方向）。
    # 因此，β是测线方向与正东方向的夹角？
    # 公式：W = 2 * D * tan(theta_half) / (cos(alpha_rad) * cos(beta_rad))
    W = 2 * D * math.tan(theta_half) / (math.cos(alpha_rad) * math.cos(beta_rad))
    return W

# 计算某点的水深（假设线性坡度，从西向东变浅）
# 设海域中心点(0,0)深度为110m，坡度1.5°，则水深变化：D = D0 - x * tan(α)，其中x为东西方向坐标（米），西深东浅。
def depth(x, y):
    # x: 东西方向坐标（米），中心为0
    # y: 南北方向坐标（米），中心为0
    return D0 - x * math.tan(math.radians(alpha))

# 目标函数：总测量长度
def total_length(beta, d):
    # beta: 测线方向（与正东方向夹角，度）
    # d: 测线间距（米）
    # 计算需要多少条测线
    # 测线方向为β，则测线间的垂直距离为 d * sin(β) ？实际上应投影到法向。
    # 由于海域是矩形，测线数量N = ceil(宽度 / (d * sin(β))) 
    # 但更准确：海域在垂直测线方向上的投影宽度为 L_proj = L_east_west * |sin(β)| + L_north_south * |cos(β)|
    # 测线数量N = ceil(L_proj / d) + 1
    beta_rad = math.radians(beta)
    L_proj = L_east_west * abs(math.sin(beta_rad)) + L_north_south * abs(math.cos(beta_rad))
    N = int(np.ceil(L_proj / d)) + 1

    # 每条测线长度：海域在测线方向上的投影长度
    L_line = L_east_west * abs(math.cos(beta_rad)) + L_north_south * abs(math.sin(beta_rad))

    total_len = N * L_line
    return total_len, N

# 检查覆盖和重叠率约束
def check_constraints(beta, d):
    # 返回惩罚项（penalty），若无约束违反则为0
    penalty = 0.0
    beta_rad = math.radians(beta)
    # 计算测线数量N
    L_proj = L_east_west * abs(math.sin(beta_rad)) + L_north_south * abs(math.cos(beta_rad))
    N = int(np.ceil(L_proj / d)) + 1

    # 生成各测线位置（垂直于测线的坐标）
    positions = np.linspace(-L_proj/2, L_proj/2, N)  # 第一条测线位置？实际上应从边界开始覆盖。

    # 对于每对相邻测线，计算重叠率
    for i in range(N-1):
        p1 = positions[i]
        p2 = positions[i+1]
        d_actual = abs(p2 - p1)  # 实际间距

        # 计算这两条测线中点的水深？近似取中心点水深？
        # 由于坡度较小，近似用中心点水深D0？
        D_mid = D0  # 近似
        W = coverage_width(D_mid, beta)

        overlap_ratio = 1 - d_actual / W
        if overlap_ratio < 0.1:
            penalty += (0.1 - overlap_ratio) * 10000  # 惩罚漏测或重叠过低
        if overlap_ratio > 0.2:
            penalty += (overlap_ratio - 0.2) * 10000  # 惩罚重叠过高

    return penalty

# 模拟退火主函数
def simulated_annealing():
    # 初始状态
    beta_current = 45.0  # 初始方向角（度）
    d_current = 1000.0   # 初始间距（米）
    current_length, N = total_length(beta_current, d_current)
    current_penalty = check_constraints(beta_current, d_current)
    current_obj = current_length + current_penalty

    T = 1000.0
    T_min = 1e-5
    alpha = 0.99
    iterations = 10000

    best_beta = beta_current
    best_d = d_current
    best_obj = current_obj

    for i in range(iterations):
        # 生成新状态
        beta_new = beta_current + np.random.normal(0, 1.0)  # 小扰动
        d_new = d_current + np.random.normal(0, 50.0)
        if d_new < 10: d_new = 10  # 防止负间距

        new_length, N = total_length(beta_new, d_new)
        new_penalty = check_constraints(beta_new, d_new)
        new_obj = new_length + new_penalty

        # Metropolis准则
        delta = new_obj - current_obj
        if delta < 0 or np.random.rand() < np.exp(-delta / T):
            beta_current, d_current = beta_new, d_new
            current_obj = new_obj
            current_length, current_penalty = new_length, new_penalty

        if current_obj < best_obj:
            best_beta, best_d = beta_current, d_current
            best_obj = current_obj

        T *= alpha
        if T < T_min:
            break

    return best_beta, best_d, best_obj

# 运行SA
best_beta, best_d, best_obj = simulated_annealing()
print(f"最优方向β: {best_beta}度")
print(f"最优间距d: {best_d}米")
print(f"目标函数值: {best_obj}")